Disorder-induced temperature-dependent transport in graphene: Puddles, impurities, 

activation, and diffusion 
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We theoretically study the transport properties of both monolayer and bilayer graphene in the 
presence of electron-hole puddles induced by charged impurities which are invariably present in the 
graphene environment. We calculate the graphene conductivity by taking into account the non- 
mean-field two-component nature of transport in the highly inhomogeneous density and potential 
landscape, where activated transport across the potential fluctuations in the puddle regimes coexists 
with regular metallic diffusive transport. The existence of puddles allows the local activation at 
low carrier densities, giving rise to an insulating temperature dependence in the conductivity of 
both monolayer and bilayer graphene systems. We also critically study the qualitative similarity 
and the quantitative difference between monolayer and bilayer graphene transport in the presence 
of puddles. Our theoretical calculation explains the non-monotonic feature of the temperature 
dependent transport, which is experimentally generically observed in low mobility graphene samples. 
We establish the 2-component nature (i.e., both activated and diffusive) of graphene transport arising 
from the existence of potential fluctuation induced inhomogeneous density puddles. The temperature 
dependence of the graphene conductivity arises from many competing mechanisms, even without 
considering any phonon effects, such as thermal excitation of carriers from the valence band to 
the conduction band, temperature dependent screening, thermal activation across the potential 
fluctuations associated with the electron-hole puddles induced by the random charged impurities 
in the environment, leading to very complex temperature dependence which depends both on the 
carrier density and the temperature range of interest. 

PACS numbers: 72.80.Vp, 72.10.-d, 73.22.Pr, 81.05. ue 
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I. INTRODUCTION 

Graphene, as a novel gapless two dimensional (2D) chi- 
ral electron-hole system, has attracted great interest in 
recent years, both experimentally and theoretically^'™. 
Its transport properties have been at the center of key 
fundamental and technological efforts with vast poten- 
tial for applications in future nanotechnologj^. For 
monolayer graphene (MLG), the fundamental interest 
arises from its unique linear chiral Dirac carrier disper- 
sion with a zero energy gap between conduction and va- 
lence band^. The bilayer graphene (BLG) is also intrigu- 
ing as its physical properties lie between MLG and 2D 
semiconductor-based electron gas (2DEG) systems which 
are gapped and non-chiral with a quadratic band dis- 
persion. Much of the early work on graphene trans- 
port focused on the density-dependent (i.e., gate voltage 
tuned)ii^"— and temperature-dependenl^i^rJ^ conductiv- 
ity in homogeneous MLG and BLG systems. The basic 
graphene transport properties, particularly at high den- 
sities far from the charge neutral Dirac point, are now 
reasonably well- understood^. 

However, unintended charged impurities, which are in- 
variably present in the graphene environment, (e.g., the 
substrate-graphene interface), lead to the formation of 
inhomogeneous electron-hole puddles in the systeroi^iii, 
which have been confirmed by experiments^^^ using the 
techniques of scanning potential and tunneling micro- 
scopies. Although MLG samples show a metallic behav- 
ior at high densities a weak "insulating" temperature- 



dependent conductivity a(T) has been measured at 
low carrier density and at the charge neutrality point 
(CNP)2i. (We define insulating/metallic temperature de- 
pendence of conductivity cr{T) as da(T)/dT being pos- 
itive/negative at fixed gate voltage.) In addition, a re- 
cent experimentii on low mobility MLG grown by chem- 
ical vapor deposition (CVD) shows a strong "insulating" 
behavior at low temperatures and a metallic feature at 
high temperatures manifesting a non-monotonic temper- 
ature dependence in the measured electrical conductiv- 
ity. In BLG samplesi^"— the strong insulating behavior 
in the temperature dependent conductivity has been ob- 
served not only near CNP but also at carrier densities 
as high as lO^^cm"^ or higher. To be more specific, in 
Rcf. 18], cr(T) in BLG increases by 20 - 40% as temper- 
ature T increases from 4 — 300 K for carrier density in 
the range 3.19 x 10^^ - 7.16 x 10^^ cm"^^ Xo understand 
this anomalous temperature dependence in cr(T), both 
MLG and BLG, it is essential to know the role of disorder 
in graphene transport. We note that phonon scattering 
(Ref. [13113) J although being weak in graphene, always 
contributes an increasing resistivity with increasing tem- 
perature and thus always leads to metallic behavior, and 
thus cannot be the mechanism for the intriguing insulat- 
ing temperature dependence often observed in graphene 
transport at lower carrier densities - in fact, at very high 
temperatures (> 300 K) graphene should always mani- 
fest metallic temperature dependence in its conductivity 
due to phonon scattering effects which we would ignore in 
the current work. Our goal here is to theoretically study 
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in a comprehensive manner the temperature dependence 
of graphene transport properties arising entirely from the 
disorder effects. 

The experimentally measured anomalous temperature 
dependent conductivity of BLG has been theoretically in- 
vestigated by applying the analytic statistical theory to 
the inhomogeneous potential fluctuation and it is found 
that the anomalous BLG (t{T) is likely to be caused 
by the electron-hole puddles induced by randomly dis- 
tributed disorder in the graphene environments^. In this 
paper, we extend this work and apply the same analytic 
statistical theory to MLG systems and explain the in- 
triguing coexistence of both metallic and insulating fea- 
tures of MLG (T(r). In the presence of large fluctuat- 
ing potentials V{r) associated with microscopic config- 
urations of Coulomb disorder in the system, the local 
Fermi level, /i(r) — Ep — V{r), would necessarily have 
large spatial fluctuations. We carry out an analytical 
theory implementing this physical idea by assuming that 
the value of the potential at any given point follows a 
Gaussian distribution, parametrized by s = Vrms (the 
root-mean-square fluctuations or the standard deviation 
in V^(r) about the average potential). This distribution 
can then be used to average the local density of states 
to obtain effective carrier densities, which can then be 
used to compute the physical quantities of interest^. 
The observed anomalous temperature dependent <t(T) is 
then understood as the competition between the thermal 
activation of carrier density and temperature-dependent 
screening effects. Our theory explains the suppression of 
the insulating behavior in higher mobility samples with 
lower disorder, which is consistent with experimental ob- 
servations. We also provide the similarity and the quanti- 
tative difference between monolayer and bilayer graphene 
transport in the presence of puddles. 

The motivation of our theory comes from the observa- 
tion that the electron-hole puddles, which dominate the 
low-density graphene landscape, allow for a 2-component 
semiclassical transport behavior, where the usual metal- 
lic diffusive carrier transport is accompanied by transport 
by activated carriers which have been locally thermally 
excited above the potential fluctuations imposed by the 
static disorder. This naturally allows for both insulating 
and metallic transport behavior occurring preferentially 
respectively at lower and higher carrier densities since 
the puddles disappear with increasing carrier density due 
to screening. At zero temperature (where no activation 
is allowed) or at very high carrier density (where pud- 
dles arc suppressed), only diffusive transport is possible. 
But at any finite temperatures and at not too high den- 
sities, there would always be a 2-component transport 
with both activated and diffusive carriers contributing to 
conductivity. Our theory develops this idea into a con- 
crete description. We emphasize that our theory explic- 
itly takes into account the inhomogeneous nature of the 
graphene landscape and is non-mean-field as a matter of 
principle. 

This paper is organized as follows. In Sec. |lTl 



we introduce the analytical statistical theory to de- 
scribe random electronic potential fluctuations created 
by charged impurities in the environment. We also cal- 
culate the modified density of states and the correspond- 
ing temperature-dependent effective carrier density in 
monolayer graphene. Then, in Sec. IIIIl we describe the 
calculations and the main features of the temperature- 
dependent conductivity of MLG in the presence of den- 
sity inhomogeneity. In Sec. lIVI and fVl we elaborate and 
extend our earlier results for the interplay between den- 
sity inhomogeneity and temperature in bilayer graphene 
(BLG) transport. We further discuss the connection of 
our theory to earlier theories in Sec. |VIl We discuss the 
similarities and quantitative differences among the effects 
of inhomogeneity (i.e., the puddles) on MLG and BLG 
transport and summarize our results in Sec. I VIII In Ap- 
pendix 1^ we discuss a microscopic theory to calculate 
the effects of potential fluctuation on graphene systems, 
providing a self-consistent formulation of graphene den- 
sity of states in the presence of random charged impu- 
rities near graphene/substrate interface, showing in the 
process that this microscopically calculated density of 
states agrees well with the model density of states ob- 
tained from the Gaussian fluctuations. 



II. TEMPERATURE DEPENDENT CARRIER 
DENSITY FOR INHOMOGENEOUS MLG 

It is well known that MLG breaks up into an inho- 
mogeneous landscape of electron-hole puddles, especially 
around the charge neutral point (CNP) ^^'^^i^^ . Below we 
derive an analytic statistical theory taking account of the 
effects of inhomogeneous density in monolayer graphene 
(MLG) to explain the nonmonotonic temperature depen- 
dent transport observed in MLG^^. We start by assum- 
ing that charged impurities, located in the substrate or 
near the graphene, create a local electrostatic potential, 
which fluctuates randomly about its average value across 
the surface of the graphene sheet. The potential fluc- 
tuations themselves are assumed to be described by a 
statistical distribution function P{V) where V — V{r) is 
the fluctuating potential energy at the point r = (x, y) 
in the 2D MLG plane. We approximate the probability 
P{V)dV of finding the local electronic potential energy 
within a range dV about ^ to be a Gaussian form, i.e., 

P{V)^^L=eM^V'/2s% (I) 
V 27rs^ 

where s is the standard deviation (or equivalently, the 
strength of the potential fiuctuation), which is used as 
an adjustable parameter to tune the tail-width^^. In the 
Appendix, we provide a microscopic approach to self- 
consistently solve the strength of potential fluctuations in 
the presence of charged impurities. Due to the electron- 
hole symmetry in the problem, we only provide the for- 
malism and equations for electron like carriers and the 
hole part can be obtained simply by changing E to —E. 
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The potential fluctuations given by Eq. ([T]) affect the 
overall electronic density of states (DOS) in MLG. In our 
model we do not assume that the size of the puddles to be 
identical, but we take the puddle sizes to be completely 
random controlled by the distribution function given in 
Eq. [T] We emphasize that our assumption of a Gaussian 
distribution for the potential fluctuations, equivalently 
implying a Gaussian distribution for the density fluctu- 
ations associated with the puddles, is known to be an 
excellent quantitative approximation to the actual nu- 
merically calculated puddle structures in graphen o^'^^ . 
The characteristics of the puddles are determined by both 
the sign and the magnitude of — Ep, i.e., a negative 
(positive) V — Ep indicates an electron (hole) region. A 
different approach utilizing equal size puddles with a cer- 
tain potential V has been used to calculate transport co- 
efficients using a numerical transfer matrix technique^. 
Then in the presence of electron-hole puddles the den- 
sity of states is increased by the allowed electron region 
fraction and given b y27i29,30 

rE 9s9v{E-V) 



De{E)=Jl 

= i?i[|erfc(- 



2TT{nVp)^ 

E , s 

-) + 



■P{V)dV 



; exp(- 



2 ^/2s' ' V27r 

where erfc(x) is the complementary error function. 



erfc(a;) = 



dt. 



(2) 



(3) 



gs9v 



^"^^^ = 2.(^.^)2 
locity, Qs =2 and 



where vp is the graphene Fermi ve- 
= 2 are the spin and valley degenera- 



cies, respectively. We have Di ^ 1.5 x 10 cm~ /meV 
with the Fermi velocity vp = 10^ m/s. Note that the 
tail of the DOS is determined by the potential fluctu- 
ation strength s. For the case s — 0, the system be- 
comes homogeneous and Df.{E) = DiE. In this case 
there is no carrier density at Dirac point {E = 0) at 
zero temperature. It is apparent that in the presence of 
potential fluctuations, the De{E) starts at finite value 
at E = and approaches DiE in high energy limit. 
For high-energy limit, the carrier is essentially free since 
nearly every point of the system is accessible. In Fig. 
[U we show the normalized density of states as a func- 
tion of energy for both electrons and holes in MLG. We 
mention that the self-consistent microscopic theory gives 
the same structure for the density of states of graphene 
systems (see Appendix E| . 

Since monolayer graphene is a semi-metal or zero-gap 
semiconductor, the electron density at finite tempera- 
tures increases due to the direct thermal excitation from 
valence band to conduction band, which is one of the im- 
portant sources of temperature dependent transport at 
low carrier densities. Therefore, we first consider the 
temperature dependence of thermally excited electron 
density. The total electron density is given by 

dE 

' ^^^^ \piE-f.)+i ^ (4) 



D{E)/{D,s) 
2.0r 




FIG. 1: (Color online). Normalized density of states for 
both electron and hole in MLG. The solid and dashed lines 
are for the DOS in inhomogeneous and homogeneous systems, 
respectively. The electron (hole) band tail locates at i5 < 
(-E > 0), which gives rise to electron (hole) puddles ai E < 
and E >0. 



where f3 = l/fcsT and ^ is the chemical potential. At 
T — 0, fj. becomes the Fermi energy /i(T = 0) = Ep. 



A. ne{T) of MLG at CNP (Ep = 0) 

When the Fermi energy is zero (or at CNP) all elec- 
trons are located in the band tail at T = and the elec- 
tron and hole densities in the band tail are given by 



no = n,{Ep = 0) = fih^Ep = 0) = Di'- 



(5) 



Note that the electron (or hole) density in the band tails 
increases quadratically with the standard deviation s. At 
finite temperatures the behavior of ne{T) at CNP be- 



ne (T) = no 



(6) 



The leading order temperature dependence in ne{T) is 
quadratic. For homogeneous MLG (s = 0) with the 
linear-in-energy behavior of the DOS, the electron den- 
sity is given by ne{T) = — Tn~^B'^^- particular, in 
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the ballistic regime the number of propagating chan- 
nels increases due to the thermal smearing of the Fermi 
surface, which leads to the observation of an insulating 
behavior in <j{T) at CNP for high mobility suspended 
graphene samples^"—. The presence of the band tail 
does not change the quadratic temperature dependence 
in the thermal excitation when the system is at the charge 
neutral point {Ep = 0). But the inhomogeneous MLG 
has no electrons in the band tails. In Fig. UK a) we show 
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FIG. 2: (Color online) (a) The electron density of MLG at 
CNP as a function of temperature for different s. At T = 
the density is given by no = Dis^ /A. (b) The temperature 
dependent electron density of MLG at finite Ef for different 
s. For s/Ef i= the leading order behavior is quadratic, (c) 
Total electron densities (solid lines) and hole densities (dashed 
lines) of MLG as a function of Ef for two different s = 30 
meV and 70 meV. The densities at the band tails are given 
by n,(EF = 0) = nn(EF = 0) = DisV4. 



the temperature dependent electron density at CNP for 
different values of standard deviation s. 



,(r) of MLG at finite doping {Ef > 0) 



In the case of finite doping (or gate voltage), i.e., Ep ^ 
0, the electron density of the homogeneous MLG (i.e., 



s = 0) is given by 

= -A 

where F\(x) — 



EdE 



exp(/3(£;- /io)) 



(7) 



/?2 

t dt 



, and /io is the chem- 

1 + exp(i — x) 

ical potential of homogeneous MLG and is determined 
by the conservation of the total electron density. Then 
the chemical potential is given by the following relation, 

E^ ^2 

— ^ — = F\{^^ixq) — Fi{—/3^o). Using the asymptotic 
formsii of the function Fi{x) for x ^ 1 and a; ^ 1, i.e., 

^2 ^2 



F^ix) 



Fiix) 



' 12 

6 



a;ln2 



X 

T 



for \x\ < 1 



e{x) + xln(l + e-l^l) for |a-| > 1, 



(8) 

we have the asymptotic formula for the chemical poten- 
tial in both low- and high-temperature limits for homo- 
geneous MLG 



^ioiT) ~ Ep 



7r2 T ^„ 
1 ( — 2 



for T<s:Tf 



Ep Tp 
41n2^ 



(9) 



for T-»Tp. 



Then the corresponding asymptotic formula of the elec- 
tron density (Eq. [7|) are given by 



"Oe(T) 



DiE 



1 + ^^) forr«T. 



, , DiEWt^ 

noe(T) ^ ^ for T » Tp 



(10) 



Since the direct thermal excitation is suppressed due to 
the finite Fermi energy, the excited electron density at 
low temperatures (T <^ Tp) increases quartically rather 
than quadratically. But at high temperatures (T 3> Tp), 
the total electron density becomes a quadratic function 
of temperature as shown for an undoped MLG. 

Next, we derive the temperature dependence of ther- 
mally excited electron density in the presence of electron- 
hole puddles (s ^ 0) at finite doping {Ep ^ 0). At zero 
temperature the electron density for the inhomogeneous 
MLG can be written as: 



ne(0) 



nh{0) 



£l^[{l + S^)eric{-^J + 



DiEl 



\/2S 

[(l + S2)erfc(-i^) 



(11) 



where s — s/ Ep. The presence of electron- hole puddles 
does not induce any additional charge in the MLG system 
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and the net carrier density n = — rih should be con- 
served. Then, the finite temperature chemical potential 
/i(T) changes as a function of both temperature and the 
strength of potential fluctuation s, and it should satisfy 
the following relation: 



JZDe{E)dE ~JgD,iE)dE 

De{E)dE 



POO 

J —c 



1 + exp{/3{E ~ /i)) 



fOO 

J — c 



Dh{E)dE 



f exp(/3(Ai - E)) ' 
(12) 

where De{E) is the electronic density of states given by 
Eq. [5] and Dh{E) = De{—E) is the density of states for 
holes. The asymptotic analytical formula of the chemical 
potential for inhomogeneous MLG is obtained as: 



/i(T) ~ Ef 



1 



7r2 T 



for T <^Tf 



(13) 



^i{T) ~ EfB{s, t) for T:^Tf 
where functions A(s) and B{s) are given as follows: 



A{S) 



^erf 



1 



V2s 



B 



(S,t)=( 



(14) 

where t = T/Tf and erf(x) — e^* dt is the error 

function. 

Combining Eqs. [2ll4land[T3l we obtain the asymptotic 
analytical formula of the electron density for inhomoge- 
neous MLG at low- and high-temperature limits as: 



ne{T)c^n,{0)+DiEl^^{l-A{S)) for T^Tf 



15) 



,{T) noe{T) + n,{0) - 



D^Ej 



for T > Tp- 



In the low temperature limit (T <C Tp), the leading order 
term for the electron density has the same quadratic be- 
havior as in undoped homogeneous MLG {Ef = 0), but 
the coefficient is strongly suppressed by fluctuation for 
the case of s < Ef, i.e., the high carrier density sample. 
While in the case of s > Ef, i.e., the low carrier den- 
sity sample, the existence of electron- hole puddles gives 
rise to a notable quadratic behavior for electron density 
ne(T) [see Fig. Hb)]. 



III. CONDUCTIVITY OF INHOMOGENEOUS 
MLG 

In this section, we calculate the finite temperature con- 
ductivity for inhomogeneous MLG with the temperature- 
dependent effective carrier density derived above. The 



existence of electron-hole puddles allows that the current 
flows through "percolation channels" and the transport 
properties of the inhomogeneous MLG system can be de- 
rived using the self-consistent effective medium theory of 
conductance in composite mixturesSi, where the num- 
ber of electrons per puddle is not an important issue for 
our theory. The percolation assumption is valid as long 
as the potential fluctuation is larger than the thermal 
energy of the carriers. Otherwise transport due to dis- 
order scattering dominates. We emphasize that in our 
formalism the crossover from the percolation transport 
to ordinary scattering-dominated diffusive transport is 
guaranteed as the temperature is increased since we are 
explicitly taking into account both diffusive transport of 
free carriers and activated transport of the classically- 
localized carriers in our theory. The only effects we ne- 
glect are quantum tunneling through the potential bar- 
riers and quantum interference since ours is a semiclassi- 
cal theory. We also do not consider Klein tunneling ex- 
plicitly in this paper because the Klein tunneling occurs 
at zero temperature for normal incident carriers at the 
electron-hole puddle boundary. We also apply the Boltz- 
mann transport theory, where we include the scattering 
mechanism with screened Coulomb impurities and short- 
range disorder—. Note that the application of Boltzmann 
transport theory is justifiable because the quantum in- 
terference effects are not experimentally observed in the 
temperature regime of interest to us in this work. It is 
conceivable that quantum interference and localization 
play some roles in graphene transport at very low tem- 
peratures, which is beyond the scope of this paper. We 
also neglect all phonon effects in this work since electron- 
phonon coupling is weak in graphene. Phonon effects are 
relevant at high temperatures (> 100 K) and have been 
considered in the literature^^'^''. 

At CNP {Ep = 0) electrons and holes are equally oc- 
cupied. As the Fermi energy increases, more electrons 
occupy increasingly larger proportion of space. As the 
Fermi energy increases to Ep 3> s, nearly all space is 
populated by the electrons [see Fig.[2jc)] and the conduc- 
tivity of the system approaches the characteristic of the 
homogeneous material. Thus, there is a possible coexis- 
tence of metallic and thermally- activated transport in the 
presence of electron- hole puddles. When electron puddles 
occupy more space than hole puddles, most electrons fol- 
low the continuous metallic paths extended throughout 
the system, but it is possible at finite temperatures that 
the thermally activated transport of electrons persists 
above the hole puddles. On the other hand, holes in hole 
puddles propagate freely, but when they meet electron 
puddles, activated holes conduct over the electron pud- 
dles. Carrier transport in each puddle is characterized by 
propagation of weak scattering transport theory'^'*. The 
activated carrier transport of prohibited regions, where 
the local potential energy V is less (greater) than Fermi 
energy for electrons (holes) , is proportional to the Fermi 
factor. If (Te and ah are the average conductivity of elec- 
tron and hole puddles, respectively, then the activated 
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conductivities are given by 



<^k\v) = ahexp[(3{V-EF)], 



(16a) 
(16b) 



where the density and temperature dependent average 
conductivities (ce and ah) are given within the Boltz- 
mann transport theory^ by Ce c>c Ueir) and ah (x nh{T), 
where Ue and Uh are average electron and hole densities, 
respectively, and (r) is the average transport relaxation 
time which includes the thermal smearing effects and de- 
pends explicitly on the scattering mechanism^ and it is 
given by, 

jdeD,{e)r{e){-df/de) 
J deDe{e)f{e) 

where r(e) and / = 1/(1 + e^^'^~^^) are, respectively, the 
energy-dependent transport scattering time and the finite 
temperature Fermi distribution function. Because the 
density inhomogeneity effects already been considered in 
the variation of effective carrier density, we use the DOS 
of homogeneous MLG -De(e) = -Die in Eq. [17] to avoid 
double counting. And r(e) is given by 



^nud^s J j^;^\{Vpk,pk')\'^gi9kk') 

X [1 - cos 6lkk'] S{epk' - epk) (18) 



T(epk) 



where epk = phvFk is the carrier energy for the pseu- 
dospin state "p" and k is the 2D wave vector, (T^k.pk') 
is the matrix element of the impurity disorder poten- 
tial in the system environment, ^kk' is the scattering 
angle between in- and out- wave vectors k and k', 
5(^kk') = [1 -|- cos^kk'] /2 is a wave function form fac- 
tor associated with the chiral nature of MLG (and is 
determined by its band structure). Udis is the appropri- 
ate 2D areal concentration of the impurity centers giv- 
ing rise to the random disorder potential^^. We con- 
sider two different kinds of disorder scattering mecha- 
nisms: (i) randomly distributed screened Coulomb dis- 
order for which ndis|(V^k,pk')P = ni\vi{q)/e{q)\'^ , where 
Vi{q) = 27re^/(Kg) is the Fourier transform of the 2D 
Coulomb potential in an effective background lattice di- 
electric constant k and e{q) = e{q, T) is the 2D finite tem- 
perature static RPA dielectric function^^ (Note that we 
use Hi to denote the charged impurity density); (ii) short- 
range disorder for which ndis\{Vpw,pk')\'^ = ndV^ where 
Ud is the 2D impurity density and Vq is a constant short- 
range (i.e. a (5-function in real space) potential strength. 
Note that the use of Born approximation for short-range 
disorder requires weak scattering conditionSI, which is 
verified by the disorder parameters we use in our calcu- 
lation. 

Now we denote the electron (hole) puddle as region 
'1' ('2'). In region 1 electrons are occupied more space 
than holes when Ef > 0. The fraction of the total area 
occupied by electrons with Fermi energy Ef is given by 
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FIG. 3: (Color online) <7t{T) of MLG at charge neutral point 
for different s (Eq. I22al and rie as given in Eq. [6]) . Inset shows 
the thermally activated conductivity of MLG as a function of 
temperature, where aa{T)/aa{0) = 1 + e''''''^/^erfc(/?s/\/2). 



s = 80meV 



12nmeV: 



p — P{V)dV. Then the total conductivity of region 
1 can be calculated. 



(Tl = - 
P 



Ef 



ia,+alt^)P{V)dV, 
'eric 



ah g^s^ ^ / Ef 



2p 



V2s V2 



4 ) -(19) 



At the same time the holes occupy the area with a frac- 
tion q = 1 — p and the total conductivity of region 2 
becomes 



0-2 



9 Jep 



(a\t^ + ah)P{V)dV 



2g 



^.^+/5i^.erfc ( ^ 



\^/2s V2 



The ai and a2 are distributed according to the binary 
distribution. The conductivity of binary system can be 
calculated by using the effective medium theory of con- 
ductance in mixtures^. The result for a 2D binary mix- 
ture of components with conductivity ai and CT2 is given 



(P-O) 



(cti - 0-2) + W (o-i - 0-2)2 + 



4(71 (72 

(2p-l)2 



(21) 

This result can be applied for all Fermi energy. For a 
large doping case, in which the hole puddles disappear, 
we have p = 1 and 172 = 0, then Eq. ((2T|) becomes a — 
CTi, i.e., the conductivity of electrons in the homogeneous 
system. 



7 



A. a{T) of MLG at CNP (Ef = 0) 

We first consider the conductivity at CNP {Ep — 0). 
The conductivities in each region are given by 



CTl 



0-2 



1 



1 



1 

2^ 



' " /2erfc(/3s 



(22a) 
(22b) 



where 77 = nh/n^ is the ratio of the hole density to the 
electron density. Since the electrons and holes are equally 
populated, we have p = q = \/2 and Oe — <Jh, then 
the total conductivity becomes at — ^0^02 — a\. The 
asymptotic behavior of the conductivity at low tempera- 
tures (fcsT ^ s) becomes 
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2kBT 2 {ksTf 



(23) 



The activated conductivity increases linearly with a slope 
I s as temperature increases. Typically s is 
smaller in higher mobility samples, which gives rise to 
stronger insulating behavior at low temperatures. The 
next order temperature correction to conductivity arises 
from the thermal excitation given in Eq. ([5]) which gives 
quadratic (T^) temperature corrections. Thus, in the low 
temperature limit the total conductivity at the CNP is 
given by: 



ot = fT(0) 



1 



2 ksT 



At high temperatures {ksT 3> s) we have 



Ct = (Te 



ksT 2{kBTY 



(24) 



(25) 



where the temperature dependence of has been given 
in Eq. The total conductivity due to the activation 
behavior approaches a limiting value and all temperature 
dependence comes from the thermal excitation through 
the change of the effective carrier density in the presence 
of the inhomogeneity given in Eq. ([5]). Thus at very high 
temperatures (T ^ s/ks) the MLG conductivity at the 
charge neutral point increases quadratically regardless of 
the sample quality. In Fig.[3]the temperature dependent 
conductivity has been calculated at charge neutral point, 
where the temperature dependent scattering mechanism 
can be neglected. In Ref. [131 ' about 60% increase of con- 
ductivity is observed as the temperature increases from 4 
K to 300 K. We estimate the potential fluctuation param- 
eter s ^ 80 meV for this sample based on our theoretical 
analysis as compared with the data. 



B. cr(r) of MLG at finite doping {Ef > 0) 

At finite doping [Ep > 0) the temperature dependent 
conductivities are very complex because three energies 



{Ep^ s, and ksT) are competing among them. Espe- 
cially when fcsT ^ s, regardless of Ep, we have the 
asymptotic behavior of conductivities in region 1 and 2 
from Eqs. and (PO]) . respectively, 
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(26b) 



where s = s/Ep and t = T/Tp. The leading order 
correction is linear but the coefficient is exponentially 
suppressed by the term exp(— i?|,/2s^). This fact indi- 
cates that in the high mobility sample with small s, the 
activated conductivity is weakly temperature dependent 
except at low density regimes, i.e. Ep < s. Since the 
density increase by thermal excitation is also suppressed 
exponentially by the same factor [see Eq. (|15l) ]. the domi- 
nant temperature dependent conductivity arises from the 
scattering timei, which manifests the metallic behavior. 
On the other hand, for a low mobility sample with a large 
s, the linear temperature dependence due to thermal ac- 
tivation can be observed even at high carrier densities 
Ef > s. 

In Fig. [4] we present the total conductivities of inho- 
mogeneous MLG as a function of temperature (a) for a 
fixed Fermi energy and several s and (b) for a fixed s 
and several Fermi energies. The calculations for Fig. |4] 
are all carried out for MLG on Si02 substrate (corre- 
sponding to dielectric constant k 2.5), charged impu- 
rity density rii = 10^^ cm~^ and short-ranged disorder 
strength UciVq = 2 (eV A)^. For total conductivity, the 
thermally activated insulting behavior competes with the 
temperature-dependent screening effects, where the lat- 
ter always give the metallic behavior in conductivity for 
MLG samples. When s is small, the activated behavior is 
suppressed and the total conductivity shows the metallic 
behavior. While for large value of s, i.e., the low mobility 
sample, the thermal activation overwhelms the metallic 
temperature dependence and the system manifests in- 
sulating behavior. For s ^ Ep the situation becomes 
much complex. At low temperatures, the leading order 
of the temperature dependence is linear (the second term 
in Eq. ([26])) and the total conductivity starts at weakly 
insulating behavior. As the temperature increases, the 
screening effects begin to dominant leading to the metal- 
lic behavior. As a result, the temperature evolution of 
the conductivity becomes non-monotonic and for large 
s (or low mobility samples) the nonmonotonic behavior 
can be more pronounced as shown in experiments^^. 



IV. TEMPERATURE DEPENDENT CARRIER 
DENSITY OF INHOMOGENEOUS BLG 

In the following of this paper, we extend our previ- 
ous studyS4 on the insulating behavior in metallic bi- 
layer graphene and compare it with MLG situation. The 
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FIG. 4: (Color online). Calculated total conductivity 
(Tt(r)/(Tt(0) of MLG with the following parameters: rn — 10^^ 
cm-2 and naVo = 2 (eV A)^ (a) at{T) for Ef = 120 meV 
and for different s. (b) at{T) of MLG for s = 80 meV and 
for several Ef ~ 80, 100, 120, 160 meV, which correspond to 
the net carrier densities n = Ue — Uh — 0.9 x 10^^, 1.2 x 10^^, 
1.5 X 10^^ and 2.3 x 10^^ cm-^ 




FIG. 5: (Color online). Normalized density of states for 
both electron and hole in BLG. The solid and dashed lines 
are for the DOS in inhomogeneous and homogeneous systems, 
respectively. The electron (hole) band tail locates at _E < 
(-B > 0), which gives rise to electron (hole) puddles at £ < 
and E>Q. 



energy further increases. 

Because BLG is also a gapless semiconductor like 
MLG, the direct thermal excitation from valence band 
to conduction band at finite temperatures composes an 
important source of temperature dependent transport in 
BLG. Thus, the temperature dependence of thermally 
excited electron density is first to be considered. 



A. ne(r) of BLG at CNP {Ef = 0) 



most important difference between MLG and BLG comes 
from the fact that, in the BLG, the two layers are weakly 
coupled by interlayer tunneling, leading to an approxi- 
mately parabolic band dispersion with an effective mass 
about m ~ O.OSSme (mg corresponds to the bare electron 
mass) contrast to linear-dispersion Dirac carrier system 
for MLG. As done for MLG, we assume the electronic 
potential fluctuations in BLG system to be a Gaussian 
form given in Eq. [T] and this potential is felt equally by 
both layera^^. 

In the presence of potential fluctuations the density of 
states (DOS) for disordered BLG is given by De{E) = 

jf^DoPiV)dV = L>oerfc(-£;/V2s)/2, where Dq ^ 

gsgv'm/ {"^Trfi^) is the DOS in a homogeneous BLG system, 
where gs = 2 and g^ = 2 are the spin and valley degenera- 
cies, respectively. We have Dq = 2.8 x lO^'^ cm^^/meV 
assuming m = O.OSSmg. The DOS of hole can be calcu- 
lated from the following relation: Dh{E) = De{—E). In 
Fig. [5l the density of states of both electron and hole 
are shown for the inhomogeneous BLG system. In the 
presence of potential fluctuations, the electron and hole 
coexist for certain amount of regions near CNP and their 
DOS approach to the homogeneous case as the carrier 



With the help of Eq. lU we could get the total elec- 
tron density for BLG in the presence of electron-hole 
puddles. We first consider the situation at CNP, where 
all electrons are located in the band tail at T = 
and the electron density in the band tail is given by 
no = ne{Ep = 0) = Dqs /^/2t:—. Contrast to the 
quadratic dependence of s in MLG, the electron density 
in the band tail for BLG is linearly proportional to the 
standard deviation s. Unlike MLG, which has the ex- 
act formula for nQ(T) (i.e., Eq. ([B])), we could only find 
the asymptotic behavior of no(r) at finite temperatures 
for BLG. The low temperature {ksT/s <^ 1) behavior of 
electron density at CNP becomes 



ne{T) = no 



^2 



(27) 



Thus, the electron density increases quadratically at low 
temperature limit. For homogeneous BLG with the con- 
stant DOS the electron density at finite temperatures is 
given by ne(T) = Dohi{2)kBT, which has the universal 
slope Z?oln(2)fcB. The presence of the band tail sup- 
presses the thermal excitation of electrons and gives rise 
to the quadratic behavior. However, at high temperature 
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limit, the density increases linearly with the same slope 
approaching to the homogeneous system, i.e.. 



n{T) ^ Do 



\n{2)kBT+l ' 



8 (fci3r)2 



(28) 



In Fig. ini[a) we show the temperature dependent electron 
density at CNP for different standard deviations. Com- 
pared with the inset of Fig. [31 it is apparent that, even for 
the same strength of potential fluctuation s, the effects 
of thermal excitation of carrier density are much stronger 
in BLG than in MLG sample, which leads to more easily 
observed insulating behavior in BLG samples. 



B. ne(r) of BLG at finite doping (Ef > 0) 

In this subsection, we derived the total electron den- 
sity at finite temperatures for inhomogeneous BLG away 
from CNP. Contrary to MLG, we need to calculate the 
finite temperature chemical potential (i.e., Eqs. [HI and 
1131) . The charge conservation relation in both homoge- 
neous and inhomogeneous BLG gives the temperature 
independent chemical potential /x = Ep, allowing us to 
directly calculate the total effective electron (hole) den- 
sity. In the case of finite gate voltage, i.e., Ep 7^ 0, the 
electron density of the homogeneous BLG for s = is 
given by 



noe{T) = DoEF l+Hn 1 



(29) 



where t = T/Tp and Tp = Ep/ks- The thermal excita- 
tion is exponentially suppressed due to the Fermi func- 
tion at low temperatures (T ^ Tp). While at high tem- 
peratures (T ^ Tp) it increases linearly. In the presence 
of finite potential fluctuations (s ^ 0), the electron and 
hole density at zero temperature for the inhomogeneous 
system are given by: 



ne(0) = DoEp 



1 



-erfc 



-erfc 

















+ ' e-1/2. 






V2^ 





(30) 



where s = s/Ep and the difference of electron and hole 
density (n — He — Uh = DqEp) is independent of the 
strength of potential fluctuation s (see Fig. Efc)). At 
low temperatures (T <^ Tp) the asymptotic behavior of 
the electron density is given by 



,{T) = n,{0) + DaEp 



-l/2s 



I2V2 



(31) 



The leading order quadratic behavior of ne{T) as in un- 
dopcd BLG {Ep = 0) is strongly suppressed by potential 
fluctuation. For the situation s > Ep, the existence of 
electron-hole puddles gives rise to a notable quadratic 
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FIG. 6: (Color online) (a) The electron density of BLG at 
CNP as a function of temperature for different s. At T = 
the density is given by no — Dqs /y/2T^. (b) The temperature 
dependent electron density of BLG at finite Ep for different 
s. For s/Ef 7^ the leading order behavior is quadratic 
while at s = the density is exponentially suppressed, (c) 
Total electron densities (solid lines) and hole densities (dashed 
lines) of BLG as a function of Ep for two different s = 30 meV 
and 70 meV. The linear line represents the density difference 
n = Ue — rih ~ DqEf, which linearly depends on the Fermi 
energy. The densities at the band tails are given by Ue^Ep = 
0) = nniEp = 0) = Dqs/V2^. 



behavior [see Fig.[6l[b)]. At high temperatures (T ^ Tp) 
we flnd 



ne(T) =noe(r) + 



DqEp Tp 

(1 + 6/3^^)2 y^- 



(32) 



where the linear temperature dependence of electron den- 
sity is dominant as the homogeneous system. 
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V. CONDUCTIVITY OF INHOMOGENEOUS 
BLG 

With the help of total electron and hole density calcu- 
lated above, we will derive the temperature dependent 
conductivity for BLG in the presence of electron-hole 
puddles. We will apply both Boltzmann theory^ and 
effective medium theory'^* to interpret the intriguingly 
insulating behavior observed in BLG sampleai^ii^i^i. 

The density and temperature dependent average con- 
ductivities in BLG, denote as a,, and are given within 
the Boltzmann transport theory: 



m 

o2/ 



(33) 



where Ue and Uh are average electron and hole densities, 
respectively, (r) is the transport relaxation time for bi- 
layer graphene: 



SdeD,{e)eT{e){^dflde) 
SdeD,{e)f{e) 



(34) 



and r(e) is calculated with Eq. [T8l But for BLG sys- 
tems, one needs to use the parabolic dispersion relation 
Epk = ph^k"^ /2m for the pseudo-spin state "p" and the 
static dielectric screening function derived in Ref. [ssj . 
The wave function form factor associated with the chiral 
nature of BLG is also different from the case in MLG, 
which is given by 5(6'kk') = [1 -|- cos 20kk'] /2. To de- 
termine the average scattering time in BLG, we take 
into account the long-range charged impurity scatter- 
ing and short-range defect scattering, which has been 
established that both contribute significantly to bilayer 
graphene transport properties^. The activated conduc- 
tivities should also be included in the presence of density 
inhomogeneity in the BLG, which follow the same rela- 
tion as given for MLG : 



e 



{V) = a,eMli{EF-V)], 
{V) = aneMPiV - Ef)1 



(35a) 
(35b) 



A. air) of BLG at CNP 



When electron-hole puddles form in the BLG samples 
(denote the electron (hole) puddle as region '1' ('2')), 
the transport properties can be treated with effective 
medium theory as described in Sec. IIIII And Eqs. 
[23] for the inhomogeneous MLG also apply to the inho- 
mogeneous BLG system. We will first discuss the total 
conductivity of BLG at CNP {Ep =0). In this case, 
the electron and hole are equally occupied and the total 
conductivity at — ci (see Eq. I22al and 1^5)) . At low tem- 
perature limit (T ^ s/ks), the activated conductivities 
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FIG. 7: (Color online) otiT) of BLG at charge neutral point 
for different s (Eq. [22a] with for BLG). Inset shows the 
thermally activated conductivity in BLG as a function of tem- 
perature, where aa{T) / aa{Q) = 1 + e'^^'''' ''^eric{l3s/^/2), the 
same as for the MLG case. 



increase linearly with a slope \j2/'KkB/ s as the tempera- 
ture increases. The next order temperature correction to 
the conductivity is quadratic T^, which arises from the 
thermal activation (see Eq. Wf\ . Thus, at low tempera- 
ture limit the total conductivity at CNP is given by 



at{T) = (7(0) 



1 



2 ksT TT^ {ksT^ ^ 
s 6 



(36) 



At high temperatures (fc^T ^ s), the total conductivity 
is given by: 



irksT 2{kBTy 



(37) 



It is apparent that the activation behavior approaches 
a limiting value at high temperature limit (T s/ks) 
while the thermally activated electron density becomes 
dominant, which increases linearly with a universal slope 
ln(2) regardless of the sample quality. Thus, all temper- 
ature dependence of the total conductivity comes from 
the thermal excitation through the change of the carrier 
density given in Eq. ([28]) . In Fig. [7] we show the cal- 
culated temperature dependent conductivity at charge 
neutral point. The inset present the activated conduc- 
tivity versus the temperature. In Ref. [l8j| , at (T) at the 
CNP of the BLG sample increases almost two times as 
temperature T varies from 4 K to 300 K. Our theoretical 
analysis using a potential fiuctuation parameter s ~ 40 
meV gives reasonable agreement with the experimental 
data. 



B. a{T) of BLG at finite doping {Ef > 0) 

The temperature dependent conductivities at finite 
doping (Ep > 0) are very complex because three en- 
ergies {Ep, s, and fcsT) are competing. Regardless of 
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Ep, when ksT <^ s, we have the asymptotic behavior of 
conductivities in region 1 and 2 the same as MLG situ- 
ation, given in Eq. (|26|) . But the average electron and 
hole conductivities ((Tg and ah) are quite different from 
MLG case, which is determined by the specific band dis- 
persion relation and also the dielectric function e{q,T). 
Thus, the leading order correction to at in BLG is also 
linear, which comes from the activated conductivity, but 
the coefhcient is exponentially suppressed by the term 
exp{—Ep/2s^). In the high mobility sample with small 
s, the activated conductivity is weakly temperature de- 
pendent except around GNP, i.e. Ep < s. Since the 
density increase by thermal excitation is also suppressed 
exponentially by the same factor [see Eq. (j3T|) ] the domi- 
nant temperature dependent conductivity arises from the 
scattering mechanism^. On the other hand, in the low 
mobility sample with large value of s, the linear tem- 
perature dependence due to thermal activation can be 
observed even at high carrier densities Ep > s. 

In Fig. [5] we calculate the total conductivities (a) for a 
fixed Ep and several s and (b) for a fixed s and several 
Ep. Even for homogeneous BLG, there are two scatter- 
ing mechanism competing with each other. The short- 
range disorder in BLG contributes to a strong insulating 
transport behavior for all temperature, whereas screened 
Coulomb scattering always leads to a metallic behavior 
for T ^ Tp2^. At low temperature limit, the total con- 
ductivity at(T) decreases with increasing temperature, 
but at higher temperatures, the short-range disorder con- 
tribution becomes quite big and leads to a at (T) increas- 
ing with T. Therefore, when s is small, the scattering 
mechanism is dominant and the total conductivity mani- 
fests a non-monotonic temperature dependence (see Fig. 
[HKa)). However, for large s the activated temperature de- 
pendence behavior overwhelms the metallic temperature 
dependence, and the system shows insulating behavior 
(see Fig. [5fb)). It clearly shows that the insulating be- 
havior in BLG sample appeared at carrier densities as 
high as 10^^ cm^^ or higher. 
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FIG. 8: (Color online). Calculated total conductivity 
at{T) / at{Q) of BLG with the following parameters: Ui — 10^^ 
and mVi = 2 (eV A)^. (a) at[T) for Ep = 60 meV 
and for different a. (b) utiT) for s = 40 meV and for several 
Ep = 20, 40, 60, 80 meV, which correspond to the net carrier 
densities n = rie - = 0.55 x 10^^ 1.1 x 10^^ 1.6 x 10^^, 
and 2.2 x 10^^ cm"^. 



VI. CONNECTION TO EARLIER THEORIES 



We have demonstrated theoretically that the observed 
insulating behavior in temperature-dependent monolayer 
and bilayer graphene conductivity can be explained by 
the thermal activation between puddles. There are also 
other theories which have been elaborated to explain low 
carrier density graphene transpor t^^'^^i'^^ . In this section, 
we establish the bridge to connect our current theory 
and earlier theories on graphene transport due to the 
formation of inhomogeneous electron-hole puddles near 
the charge neutrality point. 

The key qualitative difference between our theory and 
all earlier graphene transport theories is the introduc- 
tion of the 2-component transport model where regular 
diffusive metallic carrier transport coexists with local ac- 
tivated transport due to activation across potential fluc- 



tuations in the puddles. Our theory just explicitly ac- 
counts for the inhomogeneous landscape in the system, 
which earlier theories ignored. This 2-component nature 
of graphene transport, where both metallic and insulat- 
ing behavior coexist because of the existence of puddles, 
produces the experimentally observed complex temper- 
ature dependence with the low-density behavior being 
primarily insulating-like and the high-density behavior 
being primarily metallic-like. 

Two different theories have been developed to study 
the low-density transport in graphene, where the strong 
density inhomogeneity is dominated. In Ref. jsoj Adam 
et al. qualitatively explained the plateau-like approx- 
imate nonuniversal minimum conductivity at low car- 
rier density observed in monolayer graphene samples. 
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The basic idea is to introduce an approximate pinning 
of the carrier density at n — n* « n,j at low carrier 
density limits \n\ < \ni\, where rii is an impurity den- 
sity. The constant minimum conductivity is then given 
by (Jmin (^{n = Hi) for n < Hi. This simple the- 
ory for monolayer graphene transport qualitatively ex- 
plained the existence of conductivity minimum plateau 
and the extent to which the minimum conductivity is 
not universal, which was in good agreement with the ob- 
served density-dependent conductivity over a wide range 
of charged impurity densities^'^. However, this theory did 
not take account of the highly heterogeneous structure 
near charge neutrality point and the thermally activated 
conductivity at finite temperatures, which then can not 
explain the observed non-monotonic temperature depen- 
dent transport in low mobility graphene samples"^^. 

A more elaborate Thomas-Fermi-Dirac (TFD) theory 
and an effective medium approximation (EMT) have 
been introduced in Refs. |26|] and [l3| to study the 
electrical transport properties of disordered monolayer 
graphene. The ground state-density landscape n(r) can 
be obtained within this TFD approach and the resul- 
tant electrical transport can be calculated by averaging 
over disorder realizations and the effective medium the- 
ory. This theory gives a finite minimum conductivity and 
is able to explain the crossover of the density-dependent 
conductivity from the minimum value at the Dirac point 
to its linear behavior at higher doping. Later, this TFD- 
EMT theory is also applied to calculate the conductivity 
of disordered bilayer graphene in Ref. [3^. The TFD- 
EMT technique successfully explains the graphene, both 
MLG and BLG, transport properties in the theoretically 
difficult inhoniogcneity-dominant regime near the charge 
neutral point, but this approach fails to explain the tem- 
perature dependence of the conductivity for a wide range 
of temperatures. 

In our current model discussed above, we include three 
effects, the electron-hole structure formation, the ther- 
mal activated conductivities and the temperature depen- 
dence of screening effects, to explain the temperature- 
dependent conductivity in both monolayer and bilayer 
graphene systems. The nonmonotonic temperature- 
dependent conductivity in graphene systems is then natu- 
rally understood from the competition between the ther- 
mal activation of charge carriers and the temperature- 
dependent screening effects. Our transport theory qual- 
itatively explains the observed coexisting metallic and 
insulating transport behavior in both MLG and BLG sys- 
tems. For low mobility MLG samples, the dominant role 
on graphene conductivity switches from the thermally 
activated transport of inhomogeneous electron-hole pud- 
dles to metallic temperature-dependent screening effects, 
which gives rise to a nonmonotonic behavior from the 
strong insulating behavior at low temperatures to metal- 
lic behavior at high temperatures. On the other hand, 
another nonmonotonic temperature-dependent transport 
can be observed in very high mobility bilayer graphene 
devices, i.e., from metallic behavior at low temperatures 



due to the screening effects of Coulomb scattering to in- 
sulating behavior at high temperatures due to the short- 
range disorder. The merit of our model is that it is so 
simple that we could get the asymptotic behavior at low 
and high temperature limits analytically. Moreover, it 
provides a clear physical picture of the dominant mech- 
anisms at different regimes as discussed above. 



VII. DISCUSSIONS AND CONCLUSIONS 

We first discuss the similarity and the difference be- 
tween MLG and BLG transport from the perspective of 
our transport-theory considerations. We find that both 
manifest an insulating behavior in at (T) for low mobility 
samples. We also find that both systems could exhibit a 
non-monotonic temperature dependent conductivity for 
low mobility samples. However, the physical origin for 
the non-monotonic temperature dependence is quite dif- 
ferent in the two systems: in the MLG the non-monotonic 
feature comes from the competition between thermal ac- 
tivation and the metallic screening effects, which leads to 
at (T) first increasing and then decreasing with increasing 
temperature (see Fig. Ufa)). While for BLG, the com- 
petition between short-range insulating scattering and 
metallic Coulomb screening effects leads to at{T) first 
decreasing and then increasing as temperature increases 
(see Fig. EJa)). Most important quantitative difference 
between MLG and BLG transport comes from their band 
dispersions, which leads to much weaker effects of den- 
sity inhomogeneity in MLG so that the anomalous insu- 
lating temperature dependence of a{T) is typically not 
observed in MLG away from the CNF although the gate 
voltage dependence of MLG and BLG conductivities are 
similar The linear Dirac carrier system for MLG 
leads to linear DOS, which goes to zero at CNP, but the 
parabolic band dispersion relation in BLG leads to a con- 
stant DOS. Due to the difference in the density of states 
between homogeneous MLG and BLG, the modified DOS 
in inhomogeneous MLG is increased (see Fig. [T]) rather 
than decreased in inhomogeneous BLG (see Fig. [5]). The 
dimensionless potential fluctuation strength s (= s/Ep) 
is much weaker in MLG than in BLG from simple es- 
timates: sblg/smlg ~ 32/v^ where n = n/10^°, and 
sblg ^ Smlg upto n — IQ^^ cm^^. Direct calculationsi 
show that the self-consistent values of s tend to be much 
larger in BLG than in MLG for identical impurity dis- 
order. In addition, the qualitatively different DOS leads 
to much stronger effective short-range scattering in BLG 
compared with MLG even for the same bare scattering 
strength. Thus, the insulating behavior in at{T) will 
show up at high temperatures even for relatively higher 
mobility BLG samples (i.e., small s). In contrast, only 
in very low mobility MLG samples, where s is very large, 
can the insulating behavior of temperature dependent re- 
sistivity be observed^iil. No simple picture would apply 
to a gapped (A^) BLG system, since four distinct energy 
scales {s, Ep,kBT, and A^) will compete and the con- 
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ceivable temperature dependence depends on their rel- 
ative values^i^i^. Our assumption of BLG quadratic 
band dispersion is valid only at low (< 5 x 10^^ cm^'^) 
carrier densities, where most of the current transport ex- 
periments are carried out. At higher densities the band 
dispersion is effectively linear and the disorder effects on 
Ut{T) are weaker. 

Before concluding, we emphasize that our theory is 
physically motivated since puddles are experimental facts 
in all graphene samples. Puddles automatically imply a 
2-component nature of transport since both diffusive car- 
riers and activated carriers can, in principle, contribute 
to transport in the presence of puddles. Of course, 
the effect of puddles is much stronger at low carrier 
densities, explaining why insulating (metallic) temper- 
ature dependence is more generic at low (high) graphene 
carrier densities. We emphasize that local carrier ac- 
tivation in puddles is just one of (at least) four dif- 
ferent independent transport mechanisms contributing 
to the temperature dependent conductivity. The other 
three are temperature dependent screening (Ref. [Tl|). 
phonons (Refs. jS^H^), and Fermi surface thermal aver- 
aging (Refs. [llll33| ). Our theory presented here includes 
the three electronic mechanisms for temperature depen- 
dence: screening, Fermi surface averaging, and puddle 
activation. We leave out phonons, which have been con- 
sidered elsewhere (Ref. [2 2".' 2 3*] ) and will simply add to 
the temperature dependent resistivity. The weak phonon 
contribution to graphene resistivity makes it possible for 
the electronic mechanisms to dominate even at room tem- 
peratures, but obviously at high enough temperatures, 
the system will, except perhaps at the lowest densities 
around the CNP, manifest metallic temperature depen- 
dence with the resistivity increasing with temperature 
because of phonon scattering. Similarly, the puddle ef- 
fects dominate low densities and therefore, the insulating 
behavior will persist to very high temperatures around 
the zero-density CNP since activation across potential 
fluctuations are dominant at the CNP. It is gratifying 
to note that these are precisely the experimental obser- 
vations. We note that in general the temperature de- 
pendent conductivity of graphene could be very com- 
plex since many distinct mechanisms could in princi- 
ple contribute to the temperature dependence depend- 
ing on the carrier density, temperature range, and dis- 
order in the system. Inclusion of phonons (at high tem- 
peratures) and quantum localization (at low tempera- 
tures) effects, which are both neglected in our theory, can 
only complicate things further. What we have shown in 
this work is that the low-density conductivity near the 
CNP is preferentially dominated by density inhomogene- 
ity and thermal carrier activation effects leading to an 
insulating temperature dependence in the conductivity 
whereas the high-density conductivity, where the puddles 
are screened out, is dominated by a metallic conductiv- 
ity due to temperature-dependent screening effects. This 
general conclusion is consistent with all experimental ob- 
servations in both MLG and BLG systems to the best of 



our knowledge except at very high temperatures where 
phonon effects would eventually lead to metallic behavior 
at all densities. 

To conclude, we have investigated both MLG and BLG 
transport in the presence of electron-hole puddles within 
an analytic statistical theory. Our theory explains the 
experimentally measured insulating behavior at low tem- 
peratures and the consequent nonmonotonic behavior for 
low mobility sample o^^i^°'^^ . A reasonable quantitative 
agreement with the experimental data can be obtained 
by choosing appropriate disorder parameters in our the- 
ory (i.e. potential fluctuation and impurity strength) for 
different samples. We find that the puddle parameter 
s, defining typical potential fluctuations, to be around 
10 — 80 mcV in typical graphene samples as extracted 
by fitting our theory to existing experimental transport 
data near the charge neutrality point. These values of 
potential fluctuations characterizing the graphene charge 
neutrality point are very consistent with direct numerical 
calculations of graphene electronic structure in the pres- 
ence of quenched charged impurities^ '^"^i^^i^^ . We also re- 
late our current model to earlier theories using the picture 
of diffusive transport through disorder-induced electron- 
hole puddles. Finally, we show the similarity and the 
quantitative difference between MLG and BLG transport 
in the presence of puddles. 
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Appendix A: A self-consistent formulation of 
graphene density of states in the presence of 
inhomogeneity 

Below we provide a microscopic theory to calculate 
self-consistently the electronic density of states in the 
presence of the potential fluctuations caused by ran- 
dom charged impurities located near graphene/substrate 
interface, which has been applied to two dimensional 
semiconductor based electron gas systems^*. This self- 
consistent approach mainly addresses two problems with 
the presence of random charged impurities. One is the 
screening of the long-range Coulomb interactions be- 
tween the carriers and the charged impurities. The other 
is the real-space potential fluctuations produced by the 
random array of charged impurities. 

The motivation for this appendix is two- fold: (1) pro- 
viding a microscopic self-consistent theory of graphene 
density of states in the presence of puddles; (2) show- 
ing that our approximate physically-motivated density 
of states (Eq. ^ is an excellent approximation to the 
self-consistent density of states. 
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FIG. 9; (Color online) Standard deviation of potential fluc- 
tuation s versus the screening constant qrp (loglog plot) in 
MLG by varying the Fermi level. The dotted blue line is for 



riimp = 1.0 X 10 
solid red line is for n 
d 



zo 



1 nm and d = 100 nm. The 



imp 



0.5 X 10^"^ cm""', zo = 1 nm and 



-,12 



200 nm. The dashed green line is for Uimp = 0.5 x 10 
cm~^, Zo — 1 nm and d = 100 nm. The dotdashed black line 
is for Tiimp ~ 0.5 X 10^^ cm~^, zq — 2 nm and d = 100 nm. 



is the dielectric constant (k ~ 2.5 for graphene on Si02 
substrate). 

The screening constant shown in Eq. lAll enters Pois- 
son's equation for the potential change 0(r, z) produced 
by a charge density pext (associated with the charged 
impurities in the graphene/substrate environment). For 
a charge Ze (we use Z = 1 in the calculation) located 
at r = a/x^ + — and z — zq > (on top of the 
graphene layer), the additional Coulomb potential satis- 
fies: 



{r,z) - 2qTFg{z)<i)o{r) 



'^'nZe5(x)5{y)5{z — zq) 



(A2) 

where kq = k« = 1.0 in the vacuum (2: > 0), ko = Kins = 
3.9 in Si02 {z < 0) and k = Por graphene, 

g{z) — S{z) is the carrier density distribution normal to 
the interface and (/'o(^) = / 4>{fi z)g[z)dz = (p{r,0). 

To solve Eq. IA2I we take advantage of the cylindrical 
symmetry to write^ : 



(r,z) 



Jo{k'r)Ak'{z)k'dk' 



(A3) 
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FIG. 10: (Color online) Calculated the density of states of 
electron DeiEp) in MLG versus electron density using the 
following parameters: the insulator thickness d — 100 nm, the 
impurity distance from the interface zq — 1 nm. The solid red 
line is for unperturbed density of states. The dashed green, 
dotted blue and dotdashed black lines are corresponding to 
nimp = 0.5, 1.0 and 2.0 x 10^^ cm~^, respectively. 



1. Monolayer graphene 

First, we apply the self-consistent consideration of ran- 
dom charged impurities on the density of states in mono- 
layer graphene. 

The simple theory of linear screening gives^: 



27re" 



DeiEp) 



(Al) 



where qtf is the Thomas-Fermi screening wavevector, 
De{Ep) is the density of states at the Fermi level and k 



The potential will satisfy Eq. IA2I if 



k^Ak-2qTFAkiO)g{z) 



2Zed{z — zq) 



dz^ 



Kq 



(A4) 



At the interface z = 0, Ak{z) must be continuous and sat- 
isfy Ky{dAk/dz) - Kins{dAk/dz) = 2qTFAk{0)K. Akiz) 
should also satisfies the boundary condition Ak{z) — >■ 
as z — !■ 00. In addition, the impurity potential (f) 
will go to zero at the metallic contact below the Si02 
(i.e., Ak{—d) = and d is the thickness of the Si02 
layer). Such screening effects are absent in the Si02. 
After some algebra, the explicit expression of Ak{k,0) 
for insulator thickness d and the impurity distance from 
graphene/substrate interface zq is given by: 



Afc(fc,0) 



26"'=^'' sinh(dfc) 



kKinsCOsh{dk) + [knv + 2(7TFK)sinh(c?A:) 

(A5) 

For the thickness of insulator in the limit — > 00, we 

have Ak{k,Q) = — — , which has been given in 

[k + qTFjn 

the Appendix B of Ref. [45[ . The potential fluctuations 
with an array of point charges at random positions in 
the plane z — zq have a mean-square variation about the 
average potential^: 



[Ak{Q)fkdk 



(A6) 



To obtain specific results for the electronic density 
of states and the screening constant we use the simple 
Gaussian broadening approximation for the density of 
states^!. The disorder-induced potential energy fluctua- 
tions is described by P{V) = , ^ ^ exp{—V'^/2s'^) (Eq. 
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FIG. 11: (Color online) Calculated density of states of elec- 
tron De{E) of MLG versus energy E for different impurity 
configurations and carrier densities n. The solid red lines are 
for the non-interacting MLG system, (a) Calculated De{E) 
in MLG for the insulator thickness d = 100 nm, the impu- 
rity distance from the interface zq = 1 nm and carrier density 
n = 2.0 X 10^^ cm"'^. The dashed green, dotted blue and 
dotdashed black lines are corresponding to nimp = 0.5, 1.0, 
2.0 X 10^^ cm"^ respectively, (b) Calculated De{E) in MLG 
for d = 100 nm, ntmp = 0.5 x 10^^ cm"^ and n = 2.0 x 10^^ 
cm~^. The dashed green, dotted blue and dotdashed black 
lines are corresponding to zq = 1, 2, 3 nm, respectively, (c) 
Calculated De{E) in MLG for d — 100 nm, zq = 1 nm and 
n-imp ~ 0.5 X 10^^ cm~^. The dashed green, dotted blue 
and dotdashed black lines are corresponding to n = 0.5, 1.0, 
2.0 X 10^^ cm"^, respectively. 



IA6p . Then the density of states becomes 
°° 2Tr[nvF) 

(A7) 

= „.[|e*,--|-, + -|=„p(-|l,l 

where erfc(a;) is the complementary error function, s = 



Vrms, Di = — 77— -TTT, is the graphene (Fermi) veloc- 

2TT(hVFj 

ity, gs = 2 and = 2 are the spin and valley degenera- 
cies, respectively. 

By choosing the chemical potential Ep as a tuning 
parameter we have the following coupled equations: 

"iTF - —^DeiEp) (A8) 
s2 = 27m,„ipe^ J[AkiO)]'^kdk 

For fixed values of Ep, rUmp, d and zo, we get the self 
consistent results for s, qtf by solving the above two 
coupled equations. The electron density could be gotten 
from the formula: 

/CO 
D,{e)f{e)de (A9) 

where /(e) is the Fermi-Dirac distribution function. 
The electron density in the presence of disorder-induced 
electron-hole puddles has been discussed in SecUll where 
we use the potential fluctuation s as a fixed parameter. 
And here we self-consistently solve the parameter s from 
a microscopic point of view, which is in good agreement 
with the results shown in Sec. HT] The potential fiuctu- 
ation in Eq. IA6I affects the electronic density of states. 
But the fluctuations depend on the screening via Eq. IA4I 
while the screening depends on the density of states via 
Eq. lAll Therefore, we have a coupled problem which 
must be solved self-consistently. 

In Fig. [3J the standard deviation of the potential fluc- 
tuation s and the screening constant qtf are plotted for 
different values of the Fermi level. The self-consistently 
solved parameters (s, ^tf) depend on the fixed charged 
impurity density n.imp, the Si02 thickness d, the location 
of the fixed charged impurity zq, and the Fermi level Ep 
(i.e. the carrier density n). All these four effects can be 
understood from physical intuition. The reduction of the 
Si02 thickness weakens the potential fluctuations when 
the screening length is small even though there is a little 
effect for strong screening. As the charged impurities go 
away from the graphene layer the potential fluctuations 
is also reduced, while the potential fluctuations becomes 
stronger with the higher impurity density. Increasing the 
carrier density n gives rise to the stronger screening ef- 
fects, and leads to weaker potential fluctuations. 

In Fig.[TUl the density of states of monolayer graphene 
is given with the parameters of Si02 thickness d — 100 
nm and the distance of fixed charged impurities — 1 
nm for different impurity densities riimp- In Fig. 1111 
we present the electronic density of states for different 
carrier densities and impurity configurations. The self- 
consistent calculation of the density of states verifies the 
results presented in Sec|TT] as shown in Fig. [U where we 
choose the potential fiuctuation s as an adjustable pa- 
rameter. For monolayer graphene, the presence of spa- 
tially random charged impurities increases the electronic 
density of states in the whole range of energy. The corre- 
sponding hole density of states can be obtained by chang- 
ing the sign of energy Dh{E) = De{—E). 
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FIG. 12: (Color online) Standard deviation of potential fluc- 
tuation s versus the screening constant qtf (loglog plot) in 
BLG by varying the Fermi level. The dotted blue line is for 
riimp ~ 1.0 X lO'^^ cm~^, 20 = 1 nm and d = 100 nm. The 



solid red line is for Ui. 



0.5 X 10^^ cm 



zo 



1 nm and 



d = 200 nm. The dashed green line is for riimp = 0.5 x 10 
cm~^, Zo — 1 nm and d — 100 nm. The dotdashed black line 
is for riimp ~ 0.5 x 10^ 



cm , Zo — 2 nm and d = 100 nm. 
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FIG. 13: (Color online) Calculated density of states of elec- 
tron De{EF) of BLG versus electron density using the fol- 
lowing parameters: the insulator thickness d = 100 nm, the 
impurity distance from the interface zo — 1 nm. The solid red 
line is for unperturbed density of states. The dashed green, 
dotted blue and dotdashed black lines are corresponding to 
= 0.5, 1.0 and 2.0 x 10 cm , respectively. 



2. Bilayer graphene 
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FIG. 14: (Color online) Calculated density of states of elec- 
tron De{E) of BLG versus energy E for different impurity 
configuration and carrier densities n. The solid red lines are 
for the non- interacting BLG system, (a) Calculated De{E) 
in BLG for the insulator thickness d — 100 nm, the impu- 
rity distance from the interface zo — 1 nm and carrier density 
n = 2.0 X 10^^ cm~^. The dashed green, dotted blue and 
dotdashed black lines are corresponding to riimp ~ 0.5, 1.0, 
2.0 X 10^^ cm"^, respectively, (b) Calculated De{E) in BLG 
for d = 100 nm, rump = 0.5 x 10^^ cm"^ and n = 2.0 x 10^^ 
cm~^. The dashed green, dotted blue and dotdashed black 
lines are corresponding to zq = 1, 2, 3 nm, respectively, (c) 
Calculated De{E) in BLG for d — 100 nm, zo = 1 nm and 
riimp = 0.5 X lO^'^ cm~^. The dashed green, dotted blue 
and dotdashed black lines are corresponding to n = 0.5, 1.0, 
cm~^, respectively. 



In this subsection, wc provide the density of states 
in bilayer graphene in the presence of potential fluctu- 
ations. As shown for monolayer graphene, we use the 
linear screening written as^: 

27re2 

qTF = D.iEp) (AlO) 

K 

where D^iEp) is the density of states of BLG at the 
Fermi level and k is the dielectric constant and for BLG 



on Si02, K ~ 2.5. 

FoUowing the same procedure discussed for MLG, the 
disorder-induced potential fluctuation is described by the 
Gaussian form P{V) = ^ exp(— y^/2s^) and the cor- 
responding density of states can be written as (also see 
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SecHVl): 



(All) 



where erfc(a;) is the complementary error function, s 



Vrms (as given in Eq. IA6I) . Dq 



Qsgvm 



gs = 2 and 



gy = 2 are the spin and valley degeneracies, respectively. 
The main difference between MLG and BLG is in their 
density of states of non-interacting systems. The homo- 
geneous MLG system has the linear energy-dependent 
density of states while the density of states of the homo- 
geneous BLG is independent of energy, which leads to 
different Thomas-Fermi screening wavevectors. The po- 
tential fluctuation in Eq. IA6I affects the electronic den- 
sity of states in BLG. But the fluctuations depend on 
the screening via Eq. IA4l while the screening depends on 
the density of states via Eq. lAlOl Therefore, we have a 
coupled problem which must be solved self-consistently. 
In Fig. [121 the broadening parameter s and the screen- 



ing constant qtf are plotted for various Fermi levels (i.e. 
the carrier density n). As shown for MLG, the BLG 
parameters (s^qte) are also non-trivial function of the 
fixed charge density Uimp, the Si02 thickness d, the lo- 
cation of the fixed charged impurity zq, and the Fermi 
level. The different charged impurity configurations and 
carrier densities have similar effects on potential fluctu- 
ations of bilayer graphene as we discussed for monolayer 
graphene. The results for s(5tf) are also quite similar 
to that of MLG (in Fig. ^ only with small numerical 
difference. 



In Fig. [T31 the self-consistent electronic density of 
states of BLG has been calculated using Si02 thickness 
d = 100 nm and distance of fixed charged impurities 
Zq ^ I nni for different impurity densities. The higher 
impurity density changes the density of states more dra- 
matically. In Fig. I14[ we show the electronic density 
of states of BLG for different charged impurity config- 
urations and carrier densities. The existence of random 
charged impurities reduces the electronic density of states 
for > but create a band tail for i? < 0. 
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